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A numerical approach is presented which allows to cal- 
culate the equilibrium properties of Fermi systems which are 
both, strongly coupled and strongly degenerate. Based on a 
novel path integral representation of the many-particle den- 
sity operator, quantum and spin effects are included. As an 
illustration, results for the pressure and energy of an electron- 
proton plasma are presented. 



Coulomb systems continue to attract the interest of re- 
searchers in many fields, including plasmas, astrophysics, 
solids and nuclear matter, see Refs. |ffl,B for an overview. 
Stimulated by recent impressive experimental advances 
[Q], theoretical activities follow a large variety of ap- 
proaches: exactly solvable models kinetic theory 
and quantum statistics, e.g. computer simulations, 
e.g. p-wj and others. The most interesting phenom- 
ena such as metallic hydrogen, plasma phase transition, 
Fermi liquids, bound states etc., occur in situations where 
both Coulomb and quantum effects are relevant. How- 
ever, here, all present theories have practical difficul- 
ties. Quantum kinetic methods easily handle quantum 
and spin effects, but they become extremely complicated 
when correlation effects can not be treated perturba- 
tively. On the other hand, Monte Carlo (MC) and molec- 
ular dynamics (MD) simulations allow for an efficient 
treatment of strong coupling phenomena, but still have 
difficulties with incorporating quantum and spin statis- 
tics effect. 

In particular, there has been remarkable progress in ap- 
plying path integral quantum Monte Carlo (PIMC) tech- 
niques to Bose g and Fermi systems, see e.g. and 
Rcf. Jio| for an overview. However, there remains one 
major obstacle preventing efficient modelling of Fermi 
systems - the so-called sign problem, i.e. the appear- 
ance of differences of large numbers in the expressions 
for thermodynamic quantities. It is related to the anti- 
symmctrization of the fermionic N-particle density ma- 
trix which results in a sum of N\/2 positive and neg- 
ative terms in the thermodynamic functions. This led 
to a number of additional approximations, such as the 
fixed node and restricted path integral concepts. Despite 
impressive recent simulation results, see Mill and ref- 
erences therein, to overcome these assumptions remains 
one of the challenges from the fundamental (first princi- 
pal character) as well as practical point of view. 

In this work we demonstrate that these approxima- 



tions can in fact be avoided. We derive a different path 
integral representation for the density matrix p by in- 
troducing a modified set of integration variables. This 
leads to two major advantages for the thermodynamic 
functions: i) they do not contain an explict sum over 
permutations (with alternating sign) and ii) part of the 
volume and temperature dependencies of p can be elim- 
inated. As an example, we present explicit formulas for 
the equation of state and energy of a two-cmponent plas- 
mas, cf. Eqs. (||) and (^J). The numerical efficiency of our 
approach is demonstrated with results for the pressure 
and energy of a high-temperature/high-density electron- 
proton plasma over a wide range of degeneracy parame- 
ters. Our approach has important consequencies for sim- 
ulations of Fermi systems, such as dense quantum plas- 
mas, few-fermion systems in traps or quantum confined 
semiconductor structures and is the basis for devel- 
oping path integral MD simulations for Coulomb systems 
& . 

As is well known the thermodynamic properties of a 
quantum system of N particles are fully determined by 
the partition function Z and, consequently, by the density 
matrix fll31 



Z = J dqp(q,0;q,(3), 
V 



p(q,(3;q',l3') = (qf3\p\P'q'), (1) 



where p = exp{-(/3 - 0')H}, (3 = l/k B T, and 
q comprises the coordinates of all particles, q = 
{qi, q2, qi\r}- With an analytical expression for the 
density matrix given, one can use e.g. Monte Carlo meth- 
ods j5| 0,|fO| to evaluate the partition function and ther- 
modynamic quantities. However, for a quantum system, 
p is, in general, not known but can be constructed from 
its known high-temperature limit by means of a decom- 
position into 71 + 1 factors, each of which corresponds to 
the density matrix at an n + 1 times higher temperature 
& 
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where r( l+1 ) - t« = A/3 = /3/{n + 1). Eq. 
one view the N-particle state as a loop consisting of n+1 
vertices ("beads") located at intermediate coordinates 
— 1...7J. For quantum systems of bosons 
(fermions), furthermore, in Eq. (Q), the spin statistics 
theorem has been taken into account which requires to 
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perform an (anti-)symmetrization of the density matrix 
giving rise to the sum over all TV-particle permutations P 
with parity Kp, and P denotes the permutation operator. 
Obviously, for fermions (minus sign in the paranthcses), 
this sum contains N\ /2 positive and N\/2 negative terms. 
(To simplify the notation, we have not written the spin 
variables explicitly, their inclusion is straightforward). 

The unknown density matrix is efficiently com- 
puted by approximating the p's on the r.h.s. of Eq. (||) 
by the high-temperature density matrix p hT . For poten- 
tials which are bounded from below, the simplest choice 
is, e.g. fH|, 

p hT [«« , r« ; q^ , t< <+1 >] a po [gW , r« ; g( <+1 > , r^} 

x e" 



Po[g ( V (i W i+ V (i+1) ] = ;^ II e A 



(3) 



where po is the free-particle density matrix which is a 
product of single-particle density matrices, U is the po- 
tential energy, and Aa is the DeBroglie wave length cor- 
responding to k B (n + 1)T = l/A/3, A A ee 2irh 2 Af3/m. 
It is well known that for n — > oo, p' 11 " converges to the 
exact density matrix p. Notice that the factorization of 
p hT into a kinetic and potential energy term, Eq. (||), is 
an approximation also. The error made thereby is of the 
order of the variation of U on the spatial scale Aa and 
thus vanishes with n — + oo. This holds also for a repul- 
sive Coulomb potential, whereas the attractive electron- 
ion interaction has to be represented by a bounded from 
below effective pair potential |||],[lO]|J . 

After recalling the general concept of PIMC, we now 
explain our scheme. To simplify the presentation, we 
choose a binary mixture of quantum electrons and clas- 
sical protons, N e — Ni = N. The first crucial step is to 
transform the intermediate electron coordinates. To this 
end, we rewrite the partition function (|l|), now explicitly 
including the arguments of the ions, 



Q(N e ,Nj,0) 



with Q(N e ,N h /3)= dqdrd£p(q,[r],0). (4) 



Here, the notation q is retained for the ions, [r] sum- 
marizes the electron coordinates with r denoting the 
coordinates at the beginning of the fermionic loop and 
the dimensionless distances between neigh- 
boring vertices on the loop. Thus, explicitly, [r] = 
[r;r + A A £ (1) ;r + A A (£ (1) +£ (2) );...], and q,r,£® each 
are 3-/V-dimensional vectors. For the density matrix in 
Eq. (@) we have, using (§, |Q 



N 



p(q, [r],/3)=£>,(g, M,/3) 



(5) 
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where U(q, [r], 0) = XP{q) 

n 

+ —riJ2{Ui([r},f3) + Ur(q,[r},0)} , (6) 

and U l , Uf and Uf 1 denote the interaction energy be- 
tween ions, electrons at vertex "/" and electrons (ver- 
tex "I") and ions, respectively. Furthermore, <fL p = 

exp[— 7r|£p^ | 2 ] arises from the kinetic energy density ma- 
trix po of the electron with index p. Notice that, in con- 
trast to po in Eq. (||) , is independent of volume and 
temparature. Furthermore, we underline that, in con- 
trast to Eq. (Q) , Eq. (|5|) does not contain an explicit sum 
over the permutations and thus no sum of terms with 
alternating sign. Instead, the whole exchange problem is 
contained in a single exchange matrix given by 



\(r a -r b )+yZ 



(7) 



where y™ = Aa Zwb=i ■ Besides being a function of 
distance, this matrix contains all spin depencies of the 
N-electron subsystem. The calculation of det\%b a t,\ s which 
is crucial for evaluating the density matrix (|5|) is greatly 
simplified by its block structure (it contains two zero sub- 
matrices) which is explained in Fig. 1. As a result of 
the spin summation, the matrix carriers a subscript s 
denoting the number of electrons having the same spin 
projection. 

As a first example of applying our results (§, (|) to 
thermodynamic properties, we provide the result for the 
equation of state, (3p = dh\Q/dV = [a/3VdlnQ/da] a =i, 
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Here, <E>* e is the effective electron-ion pair potential, a is 
a scaling parameter for the length, a = L/Lq, (. . . | . . .) 



denotes the scalar product, and q p t, 
differences of two coordinate vectors: 



' pt — ' p 
x pt = r p 
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n, x pt = r p - q t , r pt = r p - r t + y p and 
q t + y p - The structure of Eq. (||) is obvious: 
we have separated the classical ideal gas part (first term) . 
The ideal quantum pressure in excess of the classical one 
and the correlation contributions are contained in the in- 
tegral term, where the second line results from the ionic 
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correlations (first term) and the e-e and e-i interaction 
at the first vertex (second and third terms respectively). 
The third and fourth lines are due to the further elec- 
tronic vertices and the explicit volume dependence of the 
exchange matrix, respectively. As a consequence of rep- 
resentation (||), the result for the pressure does not 
contain a sum over permutations. In fact, each of the 
sums in curly brackets in Eq. (|^) is bounded as the num- 
ber of vertices increases, n — > oo, which enables us to 
evaluate the pressure without additional approximations 
for the density matrix. 

Other thermodynamic quantities exhibit the same fa- 
vorable behavior. As a second result, we provide the 
formula for the energy, (3E = 6N/2 — /3dlixQ/dj3 



1 N f 

f3E = 6Nk B T + ^ £ dqdrdt Ps {q, 



[r},/3) x 



N, 



I3e 2 



P<t 



\q P t 



£ 



N, 



£ 



(n + l)\r pt 

/3e 2 



+ ££ 

p=i t 

N e 

E 



f3& e (\x pt \) 
n + 1 

Pe 2 (r l pt \y l pt ) 



=1 t=i 



:t (n + l)\r l pt \ — t «\„ -r ., , 



2(n+l)K t |3 



yy 

p—1 t—1 v 71 



/3 2 (x pt \y p ) d(3V e (x pt ) 



P t\ 



d\x 



P t\ 



o=l t=l 



n + 1 



n.li 



[3 Sdet|Cb 



det|C 



n,l I 
6 U 



80 



>, (9) 



where we denoted y l pt = y l p — y\. The structure of this 
result is similar to that for the pressure, Eq. (||). Again, 
each sum in the parantheses is bounded for n — > oo. 

Besides their advantageous analytical properties, ex- 
pressions (^|) and (^J) are well suited for numerical eval- 
uation using Monte Carlo techniques. For a fast gener- 
ation of a sequence of N-particle configurations (Markov 
chain) |]|| it is necessary to efficiently compute the ra- 
tio R of two exchange determinants corresponding to 
subsequent configurations of the Markov chain, R = 
d et \' l P2b 1 \new/det\'ip2b 1 \oid- The probability of accepting 
a MC configuration is proportional to the absolute value 
of R while the sign of the determinants is included in 
the weight function of each configuration. Without go- 
ing into details we only mention that this can be done 
using the inverse of the exchange matrix. Moreover, the 
inverse matrix is used for an efficient computation of the 
derivatives of the exchange determinants in Eqs. (|^) and 
(^J) which is expressed in standard manner as a sum of 
N determinants involving derivatives of ip^ 1 , Eq. (^). 

To demonstrate our numerical scheme, we consider 
as an example a two-component electron-proton plasma. 
For the potential <F ie we chose a simple approximation of 
the Kelbg potential which is the high-temperature limit, 
of the exact electron-ion interaction, e.g. [HQ- Thus, by 
increasing the number of electronic vertices n, in princi- 
ple, any desired accuracy can be achieved. To simplify 



the computations, we included only the dominant con- 
tribution in the sum over s corresponding to s = N/2 
electrons having spin up and down, respectively. (The 
contribution of the other terms is small and vanishes in 
the thermodynamic limit.) 

As a first test of our approach, we consider a mixture of 
ideal electrons and protons for which the thermodynamic 
quantities are known analytically, e.g. ||. Fig. 2 shows 
our numerical results for the pressure together with the 
theoretical curve. The agreement, up to values of the de- 
generacy parameter \ = n ^ 3 as large as 5 is remarkable. 
Even with only N = 32 electrons and protons deviations 
are rather small. One clearly sees that increasing the 
number of particles (see inset of Fig. 2) improves the 
numerical results systematically and extends the range 
of applicability to larger values of x- 

Let us now turn to the case of interacting electrons 
and protons. We have performed a series of calcula- 
tions in which the classical coupling parameter T — 
(47m e /3) 1 / 3 e 2 /47reofcT was kept constant while the de- 
generacy was varied. The results for the pressure and 
energy are presented in Fig. 3. One can see that for weak 
coupling and small degeneracy parameters, x < 0-5, ex- 
change effects are small, and QMC simulations without 
exchange (open circles) are close to our results. How- 
ever, with increasing x arL d T, the deviations are growing 
rapidly. Even stronger are the discrepancies with analyti- 
cal theories which are constructed as perturbation expan- 
sions, and thus are limited to small values of x an d T. 
Most strikingly is the decrease of the energy as a function 
of x predicted by the analytical models and QMC with- 
out exchange, which is in contrast to our results which 
show an increase for all values of T. We mention that for 
all results shown, the maximum statistical error is about 
5% and can be systematically reduced by increasing the 
length of the MC run. Furthermore, the accuracy is af- 
fected by the number n of vertices. In our calculations, 
we considered always temperatures above one Rydberg, 
for which we checked that it is sufficient to use n = 6. 

In summary, we presented a new path integral repre- 
sentation (H) for the A— particle density matrix which 
does not contain an explicit summation over permuta- 
tions. As a result, we were able to compute the pressure 
and energy of a strongly correlated plasma of degener- 
ate electrons and classical protons without additional as- 
sumptions for the density matrix. Our numerical results 
demonstrate the practical feasibility of our approach and 
open the way to a large variety of applications includ- 
ing dense hydrogen, fermion systems in condensed mat- 
ter, few-fermion systems in traps or semiconductor 
nanostructures as well as to MD simulations for fermions 
& 
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FIG. 1. Schematic of the exchange matrix ||^> a &||- It con- 
tains two zero sub-matrices, s denotes the number of electrons 
with the same spin projection. 

FIG. 2. Pressure (upper figure) and energy (lower figure) of 
an ideal plasma of degenerate electrons and classical proton. 
Theory (dashed line) is compared to PIMC simulations with 
varying particle number. 

FIG. 3. Pressure (upper figure) and energy (lower figure) of 
a nonideal plasma of degenerate electrons and classical proton. 
Lines with symbols are PIMC simulations for different values 
of r (see inset in upper figure). Remaining symbols denote 
quantum MC simulations without exchange (QMCNE) and 
analytical models (RSDWK - quantum statistical model of 
Rieman et al.), data from Ref. [16]. 
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